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Can catastrophic quenching be alleviated by separating shear and a effect? 
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The small-scale magnetic helicity produced as a by-product of the large-scale dynamo is believed to play a major role in dynamo 
saturation. In a mean-field model the generation of small-scale magnetic helicity can be modelled by using the dynamical quenching 
formalism. Catastrophic quenching refers to a decrease of the saturation field strength with increasing Reynolds number. It has been 
suggested that catastrophic quenching only affects the region of non-zero helical turbulence (i.e. where the kinematic a operates) and 
that it is possible to alleviate catastrophic quenching by separating the region of strong shear from the a layer. We perform a systematic 
study of a simple axisymmetric two-layer ofi dynamo in a spherical shell for Reynolds numbers in the range 1 < Rni < 10^. In the 
framework of dynamical quenching we show that this may not be the case, suggesting that magnetic helicity fluxes would be necessary. 
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1 Introduction 

It is widely believed that the solar magnetic cycle is driven by an aVt dynamo. The region of strong 



radial shear called the tachocline (Spiegel and Zahn 1992) at the bottom of the solar convection zone 



(SCZ) is believed to be the place where strong toroidal field is formed due to stretching of the weaker but 
diffuse poloidal field. It has been inferred from helioseismology that the tachocline is confined within a 



thin layer in the overshoot region which lies below the SCZ. This led Parker (1993) to propose the idea of 



an interface dynamo, where the shear is confined to a region with a greatly reduced turbulent diffusion, 
which also is the region of production of strong toroidal field. The helical turbulence generated due to 
convection and rotation in the layer above provides the turbulent a and a large turbulent diffusivity r/f 
The dynamo cycle being thus completed, the interface dynamo operates as a surface wave propagating 
along the boundary between strong shear and convection, which is also a region with a strong gradient in 
the turbulent diffusivity. 
In order to mimic the process of dynamo saturation, traditionally authors have used an algebraic quench- 



^^2 



^^2 



ing function like ao/(l + B /-Bgq) or even ao/(l + RraB /B^^), where ao is the unquenched value of a, 
Rja is the magnetic Reynolds number, B is the mean magnetic field and i?eq is the equipartition magnetic 
field. However, it appears that conservation of magnetic helicity plays an important role in the process of 
saturation. At large R^i, the magnetic helicity (/A • B dV) is fairly well conserved by the dynamo pro- 
ducing equal amounts of helicity in small and large scales, respectively. If the small-scale magnetic helicity 



is not able to escape out of the system, then the turbulent a effect is markedly reduced ( Pouquet et al 



1976 ) . This leads to the magnetic energy of the dynamo to be quenched such that the saturation value 
varies as R^^. Since astrophysical objects have large Rm {Rm ~ 10^ for the Sun), this strong dependence 
of the saturation field strength Ssat on i?m is referred to as catastrophic quenching. 

Even though the helicity constraint in direct numerical simulations (DNS) of dynamos with strong 
shear has been clearly identified, the results can be matched with mean-field models having a weaker 
algebraic quenching of a and turbulent diffusivity than a^ dynamos (Brandenburg et al. 2001). However, 



the empirically determined coefficients would depend on circumstances and are therefore not universal. 
Such a model would not obey magnetic helicity evolution and is therefore untenable on theoretical grounds. 
The interface dynamo model has been invoked several times as a way of getting around the persistent 



problem of catastrophic quenching (see Charbonneau„2005^ for a review). The belief is that the quenching 
function remains close to unity in the region of finite a since the toroidal field is expected to be weak 
there. However, to our knowledge, this has never been verified in a consistent manner for a range of 
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magnetic Reynolds numbers. In this paper we perform a series of calculations with mean-field aQ models, 
in spherical geometry, considering both algebraic and dynamical quenching formulations, for magnetic 
Reynolds numbers in the range 1 < i?m < 2 x 10^. An important feature of these models is that the 
region of strong narrow shear is separated from the region of helical turbulence as proposed in the Parker's 
interface model. 

In §2 we discuss the features of the aQ model used, and the formulation of dynamical a quenching. 
The results are highlighted in §3 and conclusions are drawn in §4. A part of the calculations presented 



in this paper will be discussed in a more detailed paper (Chatterjee et al. 2010). In this paper we focus 



specifically on the reality of the catastrophic quenching in dynamos with a and ^ effects operating in two 
widely separated layers. 



2 Nonlinear aft Dynamo 



2.1 The underlying mean-field model 

Our dynamo model consists of the induction equations for the toroidal component of the mean poloidal 
field potential, ^</,(r, 9), and the mean toroidal magnetic field, B^{r, 6), written in spherical geometry under 
the assumption of axisymmetry (d/d(j) = 0) . In some of the cases an additional evolution equation will be 
solved for the a effect (described in the next se ction). We have used a modified version of the publicly 
available solar dynamo code ^urycy described in Chatterjee et al. (2004) to perform these calculations. 
In this paper we have used a smoothed step profile for rj given by 



rj{r) = r]r + -rjt 



1 + erf 



de 



where re = 0.73Rq and de = 0.025/?©. We define the magnetic Reynolds number as Rm = Vt/Vr- In order 
to facilitate comparison with earlier work in Cartesian geometry (see, e.g., Brandenburg et al.||2009 >, it is 
convenient to define an effective minimal wavenumber, ki. Somewhat arbitrarily, we use ki = 2/Rq, which 
corresponds to a harmonic wave with 2 nodes spanning the full latitudinal extent between both poles. We 
also define the wavenumber of the energy-carrying eddies, fcf , corresponding to the inverse pressure scale 
height near the base of the convection zone. For all our calculations we have taken fcf = 7ki. Using the 
estimate r]t = u„ns/'ih we can express the equipartition field strength with respect to the turbulent kinetic 
energy as 



B, 



eq 



{A7rp)^/\rras = i^7rp)^/^3rjtkf. 



For algebraic quenching there is no magnetic a effect and we just have the kinetic a effect, ok, which we 
assume to be of the form 



«k(?") = -aocosl 



1 + erf 



r — r„ 



da 



1 + 5ai?V^eq) , 



(2) 



where the value of oq may be computed using the first order smoothing approximation (FOSA) as being 
equal to efr/t/cf (Blackman and Brandenburg 2002). Here, the prefactor ef is usually of order 0.1 or less 
since {u ■ a;)rms < ''^rms'^rms- The case ef = 1 means the flow is maximally helical. The term g^ is a 
non-dimensional coefficient equal to 1 or i^m depending on the assumed form of algebraic quenching 
in the models. Even though the helical turbulence pervades almost the entire convection zone, we take 
Tq = 0.771^0 and da = O.O15i?0 so that we can have a large separation between the shear layer and the 
layer where turbulence is important. Consequently we consider a differential rotation profile like that in 



^The code Surya and its manual can be obtained by writing an email to arnab@physics.iisc.ernet.in 
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Figure 1. Profiles of negative radial shear dQ/dr (dashed), a (solid) and rj (dashed-dotted) as a function of fractional solar radius 
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(3) 



where Qq = 14nHz, r^ = O.68i?0 and dw = O.OlSi^Q. The radial profiles of r/t, a and dQ/dr are plotted 
as a function of fractional radius t/Rq in Fig. [I] The region of strong radial shear is thus separated from 
the region of helical turbulence and the diffusivity has a strong gradient at a radius lying between the 
two layers. It may be noted that, in order to have supercritical dynamo action in a model with r] having 
the same radial profile as a, we must set ef ^ 1. If the strong gradient of r] lies between the two source 
regions, then we can work with ef < 1. Also the time period Tcyi of the oscillatory dynamo remains a 
reasonably small fraction of the turbulent diffusion time idifi- We can justify the profiles of r/ and a on the 
grounds of having dynamo action for a reasonable range of parameters even while avoiding any significant 
overlap between the two source regions. The formulation of the equation for the evolution of the a effect 
for dynamical quenching is described in §2.2. 



2.2 Dynamical a. quenching 



It was first shown by [Pouquet et al. ( 1976 ) that the turbulent a effect is modified due to the generation 
of small-scale magnetic helicity such that the total a effect is given by. 



Q = Ok + om 



[UJ ■ u 



p-'j ■ b) 



(4) 



where uj, u, j, b denote the fiuctuating components of vorticity, velocity, current density, and magnetic 
field in the plasma, respectively. 

For this type of quenching we use the same radial and latitudinal profiles of ax as given in Eq. ([2]), 
but without the algebraic quenching factor in the denominator, i.e. we put Qa = 0. The second term in 
Eq. (HI) is sometimes referred to as the magnetic a effect or om- It is possible to write an equation for the 
evolution of om from the equation for the evolution of the small-scale magnetic helicity density hf = a ■ b 



using the relation (Brandenburg et al. 2009), 



"M 



Vtkf 



hi. 



(5) 
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The equation for a • 5 is in principle gauge-dependent. However, under the assumption of scale separation, 
i.e. when the correlation length of the turbulence is small compared to the system size, one can define a 
magnetic helicity density of small-scale fields in a gauge-independent manner as the density of linkages 
Subramanian and Brandenburg ( 2006 ) . Using Eq. ([s]) , this leads to an evolution equation for au , 



dau 



-2vtk 



SB 

-^eq 



+ 



JXrn 



V-F„, 



(6) 



where £ and B are the mean electromotive force and the mean magnetic field. The flux F^ consists of 



individual contributions, e.g., advection due to the mean flow, the Vishniac-Cho flux (Vishniac and Cho 



2001, Subramanian and Brandenburg 2004), diffusive fluxes, triple correlation terms, etc. The effect of 



some of these individual fluxes in a spherical geometry and with both radial and latitudinal shear will be 



investigated in a paper under preparation (Guerrero et al. 2010), but for all the simulations presented in 
this paper we have put Fa = 0. 
We solve the equations for A^{r, 9), B^{r, 9) and om in a domain confined by < < vr and O.55i?0 < 



r < Rq. The boundary conditions for A^j, are given by a potential field condition at the surface (Dikpati 



and Choudhuri 1994) and A^ = at the poles. At the bottom we use the perfect conductor boundary 
condition d{rBQ)/dr = d{rB^)/dr = 0. Also B^ = on all other boundaries. We have checked that the 
results are not very sensitive to different boundary conditions at the bottom boundary mainly because the 
boundary is far removed from the dynamo region. Since F^ = 0, no derivatives of ayi need to be evaluated, 
so no boundary conditions need to be specified for qm and its evolution equation is just an initial value 
problem. We start with an initial dipolar solution where B^ is antisymmetric about the equator. 



3 Results 



To study the Rm dependence of i?sat in our model we keep all the dynamo parameters the same for all 
the runs except rj^. which we change from 2 x 10^ cm^ s~^ to 2 x 10^*^ cm^ s~^ while keeping r/t fixed 
at 4 X lO^'^cm^s"^. It may also be noted that the time period of the dynamo models (Tcyi) is fairly 
independent of the magnetic Reynolds number. 

To be able to correctly compare the dynamo models for different i?m; we have calculated the critical 
value of ao) denoted by Uc for each model. In the following we present results for ao = 2ac. We show in 
Fig. [2] the butterfiy diagrams for B^f, and om at a depth of O.72i20. It may be concluded from the butterfly 
diagram for ayi in Fig. ^p that the small-scale current helicity, and hence ajvi j is predominantly negative 
(positive) in the Northern (Southern) hemisphere. Let us denote the exponential decay time for om by ta. 
So, ta = Rm/vtkf = 4.55 X 10~^i?niidiflf- For iJin = 2 X 10^, the decay time ta ^ Tcy\ and so the system of 
equations is overdamped, as can be seen from the butterfly diagrams in Fig. [2^ as well as from saturation 
curve (dashed dotted line) in Fig. ^ Note that there are amplitude modulations of the magnetic field 
before it settles to a final saturation value. The nature of the saturation curves of the magnetic energy is 
thus strongly governed by the ratio of ta and Tcyi. 

The results of our calculations for different Reynolds numbers are plotted in Fig [3J The slopes in the 
kinematic phase are similar for all Rm within the error in the numerical determination of the critical Qc. 
The strong i?m dependence, which is reminiscent of catastrophic quenching in large Rm dynamos, can be 
easily discerned from the same figure, but is more clear in Fig. [4] where we see that the saturation energy 
decreases monotonically as a function of magnetic Reynolds number. For R^ = 2 x 10^, the code has to 
be run for 500 tdiff before the dynamo fields may start becoming 'strong' again for the case with ao = 2qc. 
Due to long computational times involved in this exercise we have not continued the calculation beyond 
60 tdis- Hence the determination of saturation magnetic energy may be inaccurate for R^ = 2 x 10^. In 
Fig.|4]we compare the case with dynamical quenching against the cases with a simple algebraic quenching 
of the form given in Eq. ^ with Qa = ^ and ga = Rm- We notice that for ga = Rm, the algebraically and 
dynamically quenched a effects seem to give similar dependences on R^. 

When the code is run longer, we start seeing changes in the parity after t > 40tdiflf for the dynamically 
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Figure 2. (a) B^{O.72Rq,0) and (b) am(O.72i?0, 9) as a function of diffusion time rjtkft for R^ = 2 X 10^ 
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Figure 3. Volume averaged magnetic energy in the domain scaled with the equipartition energy for Rni = 1 (dotte d line), Rm = 20 
solid) , Rm = 200 (dashed), ifm = 2 X 10"^ and i?m = 2 X 10^ (triangles) with dynamical quenching. Adapted from [Chatterjee et al.\ 

'mm- 



quenched system in contrast to the strong' algebraic quenching case, where the parity remains dipolar. 
However the magnetic energy and the dynamo period T^yi remain fairly constant even while the system 
fluctuates between symmetric and anti-symmetric parity at an irregular time interval. 

In Fig. [5] we show the meridional snapshots in the Northern hemisphere of the toroidal component 
of the magnetic field and om- It may be noted that the regions strongest in B^j, become progressively 
confined in the narrow shear layer with increasing R^ while om becomes stronger leading to decrease in 
i?sat • Even though om is predominantly negative in the Northern hemisphere there is a region of positive 
small-scale helicity generated just below the region where ok is finite so that contribution from the term 
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Figure 4. Volume averaged magnetic energy scaled with the equipartition energy in the saturation phase as a function of Rm for 
dynamical a quenching (trian gles +solid) and algebraic quenching with 3a = 1 (squares + dashed) and with g^, = Rm (cross + 
dashed-dotted) . Adapted from [Chatterjee et aZ.|([2010i). 



aB^ in the source S ■ B is small. This effect is similar to the one reported in Brandenburg et al. (2009). 
We also have not observed any evidence of chaotic behaviour in the range of magnetic Keynoids number 



20 < i?in < 2 X 10^ for twice supercritical a {a < 2ac)-, in agreement with Covas et al. (1997) 



4 Conclusions 



The calculations done in this paper indicate that it is not possible to escape catastrophic quenching due to 
accumulation of small-scale helicity in the domain by merely separating the regions of shear and a effect. 
The saturation value of magnetic energy decreases as ~ R^ for both dynamical quenching and the 'strong' 
(~ i?ni-dependent) form of algebraic quenching for the simple two-layer model. However, additionally we 
observe parity fluctuations for cases with dynamical quenching. It does not seem to us that there exists 
any chaotic behaviour in the time series of magnetic energy since the dynamo period and the saturation 
energy remains fairly constant. It may be possible that solar wind, coronal mass ejections, and Vishniac 
and Cho fluxes help in throwing out the small scale helicity from the Sun and thus alleviate catastrophic 
quenching. In this study we have not found any difference between the nature of the saturation curves for 
an afi dynamo and an q^ dynamo using the form of dynamical quenching given by Eq. ([6|). However, it is 
clear that the algebraic quenching formula must fail if we were to allow for magnetic helicity fluxes that 
would, under suitable circumstances, alleviate catastrophic quenching. 

We have been cautious about using dynamical quenching equation for dynamo numbers not very large 
compared to the critical dynamo numbers. We would expect that the magnetic field should affect all 
the turbulence coefficients including both a and rj. However for this analysis we have not included any 
quenching of r^t- This may be justified since such an effect could be mimicked by our simple two-layer model 
with a lower r/t in the region of production of strong toroidal fields and a higher r]i in the region of weaker 
poloidal fields. The effect of dynamical quenching on more realistic solar dynamo models having meridional 



circulation, Babcock-Leighton a effect, and diffusive helicity fluxes have also been studied (Chatterjee et 



al. 2010). 



Unfortunately the direct numerical simulations have not yet reached the modest Reynolds numbers used 
in this paper (~ lO'*) which are still much lower than the astrophysical dynamos. If B^^^^ really has an 
inverse dependence on i^m, then the solar dynamo should not be operating like it does! We may conclude 
that either we must have helicity fluxes out of the system or cross-equatorial diffusive fluxes inside the 
domain or that the ayi equation to be used beside the mean-field induction equations must be suitably 
modified for aVi dynamos. To verify if the equation for dynamical quenching works in the same way as 
in a^ dynamos, we need to perform systematic comparisons between DNS with shear and convection and 
mean-field modelling for ail dynamos. 
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Figure 5. Meridional snapshots of (a) B^/Beq and (b) am (color/grey-scale coded) for cases with different Rm starting from 20 (upper 
panel), 200 (middle panel) and 2x10^ (lower panel). Two concentric circles have been drawn at O.68-R0 and 0.77-Rq to denote the radial 
positions of the shear layer and the a effect. Solid and dashed lines denote poloidal field lines, corresponding to contours of positive 
(negative) rsinSA^. 
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